Decay of the ISS Orbit.

The ISS operates within the LEO environment and, given its size, at-mospheric drag results in noticeable orbit decay that requires periodic corrections via thrusters. Toexamine this in more detail, suppose that the ISS is at an altitude of 350 km and in a circular orbit.If the ballistic coefficient for the ISS is estimated as B= 165 kg/m^2, estimate the orbit decay over aperiod of six months.

In [1]:
#Rex Calabrese
import numpy as np
import numpy.linalg as la

import scipy as sp
from scipy import stats
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import Image

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

from matplotlib import rc
rc('text', usetex=False) 

fs_default = (11,4)
In [2]:
Image("images/hw8_diagram_1.png",width=600)
Out[2]:

Animation of the results!

In [192]:
from IPython.display import Video
Video('Drag_animation_2_nvidiah264.mp4',embed=True,width=600)
Out[192]:

Physical Parameters

Universal Gravatational Constant

$ G= 6.67408 \cdot 10^{-20} \frac{km^3}{kg \cdot s^2}$

Mass Earth and Jupyter

$m_{e} = 5.974 \cdot 10^{24} kg$

Radius of earth

$r_{e} = 6378 \cdot 10^{24} km$

?? grav constant earth

$\mu = 3.986 \cdot 10^{5} kg$

In [3]:
me = 5.974e24 #kg
G = 6.67259e-20 #km^3 · kg^-1 · s^-2
re = 6378.0   #km
mu = G*me # km^3 · s^-2

Satellite Parameters

$B=\frac{m}{C_{D} A}$

$\mathbf{a}_{d r a g}=-\frac{1}{2} \frac{\rho}{B}\left|\mathbf{v}_{r e l}\right| \mathbf{v}_{r e l}$

$\mathbf{v}_{r e l}=\mathbf{v}-\omega \times \mathbf{r}$

In [4]:
A = 1. * (1. / 1000)**2  #km
Cd = 2.  #unitless
m = 100 #kg

#Ballistic Coefficent
B = m / (Cd * A) #kg / km^2

print('Ballistic Coefficent: B = {0:0.3e} kg / km^2'.format(B))

Ballistic Coefficent: B = 5.000e+07 kg / km^2

Initial conditions

$\mathbf{r}_{0} = [5873.40, -658.52, 3007.49] km $

$\mathbf{v}_{0} = [-2.90, 4.09, 6.14] \ \ \frac{km}{s} $

In [5]:
r0 = [5873.40, -658.52, 3007.49] #km 
v0 = [-2.90, 4.09, 6.14] #km/s


u0 = np.hstack((r0, v0))
In [6]:
la.norm(r0)
Out[6]:
6631.400474296511

Density Modelling

$\text{for z < 100 km} : \rho = 1 \cdot 10^9 \frac{kg}{m^3} $

$\text{for z >= 100 km} : \rho = 2 \cdot 10^9 \cdot z ^ {-8.284} \frac{kg}{m^3}$

In [7]:
def rho_fn(z):
    return 0
In [9]:
# %% 

def rho_fn(z):
    if z < 100:
        return float(1e9)  
    if z >= 100: 
        return float(2e9 * z**(-8.284))
In [276]:
#Density Model: Interpolation Method
#Load data from MSISE90 model
import pandas as pd
MSISE90_pd = pd.read_excel ('density_tables/nasa_MSISE90.xlsx')  #kg/km^3
MSISE90=MSISE90_pd.values
def rho_fn(z):
    return np.interp(z, MSISE90[:,0], MSISE90[:,1]) / 1000  # / 1000 is a hack delete
In [11]:
fig00, ax00 = plt.subplots(figsize=fs_default)
for z in np.linspace(100,150,50):
    ax00.scatter(z,rho_fn(z),c='grey')
plt.show()
In [12]:
print("SANITY CHECK")
Image("images/orbitalDrag_0.png",width=500)

SANITY CHECK
Out[12]:
In [13]:
print('Model says')
print('ρ(400 km) = {0:0.3e} kg / m^3'.format(rho_fn(400) * 1000**-3))  #kg/m^3
#print('... looks reasonable')

Model says
ρ(400 km) = 5.566e-22 kg / m^3
In [14]:
'{0:0.5e}'.format(0.0000082195)
Out[14]:
'8.21950e-06'

Relative Velocity

In [10]:
omega_earth = [0., 0., 7.27e-5] #rad/s  

def v_rel_fn(r,v,omega_earth):
    return v - np.cross(omega_earth,r)

def a_drag_fn(r,r_mag,v):
    v_rel = v - np.cross(omega_earth,r)
    v_rel_mag = la.norm(v_rel)
    
    rho = rho_fn(r_mag - re)
    return   -0.5 * (rho / B) * v_rel * v_rel_mag
In [11]:
# %% State Vector Derivative Function

def deriv(u, dt):
    r = u[0:3]
    r_mag = la.norm(r)
    v = u[3:6]
    n = -mu / r_mag**3
    a_d = a_drag_fn(r,r_mag,v) 
    return [u[3],     #  dotu[0] = u[3]'
            u[4],     #  dotu[1] = u[4]'
            u[5],     #  dotu[2] = u[5]'
            u[0] * n + a_d[0],       #  dotu[3] = u[0] * n
            u[1] * n + a_d[1],       #  dotu[4] = u[1] * n
            u[2] * n + a_d[2]]       #  dotu[5] = u[2] * n

Two Body Model

In [12]:
def twobody(dt,u0):
    u = odeint(deriv, u0, dt, atol=1.0E-8) #1E-8
    r = u[:,0:3]
    v = u[:,3:6]
    return r,v
In [13]:
Y2S = 365 * 24 * 3600
D2S = 24*3600 
Y2D = 365 

#ORIGINAL: 5 years
nDays = 5 * Y2D
#SMALLER
#nDays = 1 * Y2D

ElapsedTime = nDays*D2S


TimePeriod = int(ElapsedTime / 5)

CALCULATE_HIGH_RESOLUTION = True
if CALCULATE_HIGH_RESOLUTION == True:
    
    dt_hires = np.arange(0.0, ElapsedTime, 1)  
    r1_hires = np.zeros((dt_hires.size,3))
    v1_hires = np.zeros((dt_hires.size,3))
    
    u0_i = u0
    for i in range(5): #5
        
        i_L = i * TimePeriod
        i_R = (i+1) * TimePeriod
        
        dt_i = dt_hires[i_L:i_R]
        
        r1_hires[i_L:i_R],v1_hires[i_L:i_R] = twobody(dt_i,u0_i)
        u0_i = np.hstack((r1_hires[i_R-1,:],v1_hires[i_R-1,:]))
   
    
LOAD_HIGH_RESOLUTION = False        #127 million time steps
if LOAD_HIGH_RESOLUTION == True:
    r1_hires = np.load('./results_export/full/2/r1.npy')
    v1_hires = np.load('./results_export/full/2/v1.npy')
    dt1_hires = np.load('./results_export/full/2/dt.npy')
In [263]:
r1_hires[-1]
Out[263]:
array([4394.0182087 , 5351.78300732, 1643.72185906])
In [264]:
dt_hires[-1]
Out[264]:
8639999.0
In [ ]:
In [14]:
#fix this for solving vs loading data
dt1_hires = dt_hires

Down-Sampling Data

In [15]:
# %% Chop off the inside of the high resolution trajectory data
ends = 10000
r1_ends = np.vstack((r1_hires[0:ends,:],r1_hires[-ends:,:]))
v1_ends = np.vstack((v1_hires[0:ends,:],v1_hires[-ends:,:]))
dt1_ends = np.vstack((dt1_hires[0:ends],dt1_hires[-ends:]))
NU = np.zeros(dt1_ends.size)
In [16]:
# %% Keep the tail end of hires trajectory data.

r1 = r1_hires[-100000:]  #100k
v1 = v1_hires[-100000:]  #100k
dt1 = dt1_hires[-100000:]  #100k
dt1_days = dt1 / D2S

r1_mag = la.norm(r1,axis=1)
v1_mag = la.norm(v1,axis=1)
In [17]:
# %% Downsample hires trajectory data.
sample_interval = int(D2S / 10)

r2 = r1_hires[0::sample_interval,:]
del r1_hires

v2 = v1_hires[0::sample_interval,:]
del v1_hires

dt2_days = dt1_hires[0::sample_interval] / D2S
del dt1_hires

dt2 = dt2_days / Y2D
In [303]:
# Subscript 2 data is downsampled.
# Subscript 1 data is clipped.

Calculate Orbital Elements

In [18]:
h2 = np.cross(r2,v2,axis=1)
r2_mag = np.linalg.norm(r2, axis=1)
v2_mag = np.linalg.norm(v2, axis=1)

eccvec2 = np.cross(v2, h2)/mu   -    np.divide(r2,r2_mag[:,None])
eccmag2 = la.norm(eccvec2,axis=1)


energy2 = (v2_mag**2 / 2) - (mu/r2_mag)
a2 = -mu/(2*energy2)

r2_a = a2 * (1 + eccmag2)
r2_p = a2 * (1 - eccmag2)

Plot Semi-Major axis over time

In [19]:
fig0, ax0 = plt.subplots(figsize=fs_default)
ax0.plot(dt2,a2,c='grey')

slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,a2)
line = slope*dt2+intercept
ax0.plot(dt2,line,'r--',lw=1)

ax0.set_title('Semi-Major axis :  $\Delta a$/$\Delta t$={0:0.2f} km/yr'.format(slope,intercept))  
ax0.set_xlabel('Time (yr)')
ax0.set_ylabel('Semi-Major Axis (km)')
ax0.legend(['e','linear regression'])
ax0.grid(True)
plt.show()
In [20]:
dt2[-1]
Out[20]:
4.99972602739726
In [21]:
print('semi major axis decay over {0:0.2f} years'.format(dt2[-1] ))
print('Δa = {0:0.3f} km'.format(a2[-1] - a2[0] ))

semi major axis decay over 5.00 years
Δa = -8.605 km

Plot eccentricity over time

In [22]:
fig1, ax1 = plt.subplots(figsize=fs_default)
ax1.plot(dt2,eccmag2,c='grey')

slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,eccmag2)
line = slope*dt2+intercept
ax1.plot(dt2,line,'r--',lw=1)

ax1.set_title('Eccentricity :  $\Delta e$/$\Delta t$={0:0.2e} yr$^-1$'.format(slope,intercept))  
ax1.set_xlabel('Time (yr)')
ax1.set_ylabel('Eccentricity (km/km)')
ax1.legend(['e','linear regression'])
plt.show()

Plot Perigee and Apoapsis Values over Time

In [23]:
fig3, ax3 = plt.subplots(figsize=fs_default)

c_rp = 'orange'
c_ra = 'purple'

#apoapsis
ax3.plot(dt2,r2_a - re,c=c_ra,lw=3,label='$apogee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))

#perigee
ax3.plot(dt2,r2_p  - re,c=c_rp,lw=3,label='$perigee$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p - re)
line = slope*dt2+intercept
ax3.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))


ax3.set_title('Atmospheric Drag: Perigee and Apoapsis Values')
ax3.set_ylabel('$Altitude$   (km)')
ax3.set_xlabel('Time  (years)')
ax3.legend()
ax3.grid(True)
plt.show()

Re-frame the Above Plot.

In [24]:
# %% Plot flat orbit with r_a and r_p

r2_a_shift = r2_a-r2_a[0]
r2_p_shift = r2_p-r2_p[0]

fig4, ax4 = plt.subplots(figsize=fs_default)

#apoapsis
ax4.plot(dt2,r2_a_shift,c=c_ra,lw=3,label='$r_a$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_a_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'r--',lw=1,label='$\Delta r_a$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))

#perigee
ax4.plot(dt2,r2_p_shift,c=c_rp,lw=3,label='$r_p$')
slope, intercept, r_value, p_value, std_err = stats.linregress(dt2,r2_p_shift)
line = slope*dt2+intercept
ax4.plot(dt2,line,'b--',lw=1,label='$\Delta r_p$/$\Delta t$ = {0:0.1f} km/yr'.format(slope))


ax4.set_title('Atmospheric Drag: Perigee and Apoapsis Values (arbitrary initial values)')
ax4.set_ylabel('$r$   (km)')
ax4.set_xlabel('Time  (years)')
ax4.legend()
plt.show()

The above plot shows a relatively steady perigee. We know the drag force is largest at perigee due to increased atmospheric density, similar to an instantaneous ΔV maneuver at perigee, this location remains 'pinned'.

'Manually' find r_a, r_p index locations.

In [25]:
# %% 'Micro View': FIND r_a and r_p locations, and plot over r vs t.

r1_apo = find_peaks(r1_mag)
r1_peri = find_peaks(-r1_mag)

plt.figure(figsize=fs_default)

#Plot r
plt.plot(dt1_days,r1_mag,color='k')

#Plot the found r_a,t_p locations
for i in range(r1_peri[0].size - 1): #
    j = r1_peri[0][i]
    j1 = r1_apo[0][i]
    plt.scatter(dt1_days[j],r1_mag[j], c=c_ra)
    plt.scatter(dt1_days[j1],r1_mag[j1], c=c_rp)
plt.legend(['|r|','Apogee','Perigee'],loc=1)
plt.title('Test: Find Perigee and Apoasis Values Through Peak Searching')
plt.show()

Recalculate Deceleration and Velocity Values (were not saved from inside 2body loop)

In [26]:
# %% RECALCULATE ACCEL (#Using r1, v1 ,dt1 )
a1 = np.zeros((dt1.size,3))
for i in range(dt1.size):
    a1[i,:] = a_drag_fn(r1[i,:],r1_mag[i],v1[i,:])
    
a1_mag =  la.norm(a1,axis=1)
In [321]:
a1[1]
Out[321]:
array([ 4.55131995e-19, -1.26738833e-18, -2.61900459e-18])
In [27]:
# Recalculate V_rel
v1rel = np.zeros((dt1.size,3))
for i in range(dt1.size):
    v1rel[i,:] = v_rel_fn(r1[i,:],v1[i,:],omega_earth)
    
v1rel_mag = la.norm(v1rel,axis=1)

Trying to understand where the negative deceleration is coming from.

In [28]:
dt1_yr = dt1 / Y2S
dt1_day = dt1 / D2S

fig5, ax5 = plt.subplots(3, figsize=(11,8), sharex=True)
ax5[0].plot(dt1_day, r1[:,0],c='r',label='$r_x$')
ax5[0].plot(dt1_day, r1[:,1],c='g', label='$r_y$')
ax5[0].plot(dt1_day, r1[:,2],c='b', label='$r_z$')
ax5[0].plot(dt1_day, r1_mag[:],c='k', label='|r|')

ax5[1].plot(dt1_day, v1rel[:,0],c='r', label='$v_x$')
ax5[1].plot(dt1_day, v1rel[:,1],c='g', label='$v_y$')
ax5[1].plot(dt1_day, v1rel[:,2],c='b', label='$v_z$')
ax5[1].plot(dt1_day, v1rel_mag,c='k', label='|v|')

ax5[2].plot(dt1_day, a1[:,0],c='r', label='$a_x$')
ax5[2].plot(dt1_day, a1[:,1],c='g', label='$a_y$')
ax5[2].plot(dt1_day, a1[:,2],c='b', label='$a_z$')
ax5[2].plot(dt1_day, a1_mag,c='k', label='|a|')

#r_a and r_p
for i in range(r1_apo[0].size): #
    j1 = r1_apo[0][i]
    if i ==0:
        ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp, label='A')
        ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp,label='A')    
        ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp,label='A')
    else:
        ax5[0].scatter(dt1_day[j1],r1_mag[j1], c=c_rp)
        ax5[1].scatter(dt1_day[j1],v1rel_mag[j1], c=c_rp)    
        ax5[2].scatter(dt1_day[j1],a1_mag[j1], c=c_rp)
    
    
for i in range(r1_peri[0].size): #
    j = r1_peri[0][i]
    if i ==0:
        ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra, label='P')
        ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra, label='P')   
        ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra, label='P')
    else:
        ax5[0].scatter(dt1_day[j],r1_mag[j], c=c_ra)
        ax5[1].scatter(dt1_day[j],v1rel_mag[j], c=c_ra)   
        ax5[2].scatter(dt1_day[j],a1_mag[j], c=c_ra)



for i in range(3): ax5[i].legend(loc=1)
ax5[0].set_title('Position, Relative Velocity, Deceleration (final orbit)')

for i in range(3): ax5[i].set_xlim((5*Y2D - 0.18,5*Y2D))
ax5[2].set_xlabel('Time (days)')

ax5[0].set_ylabel('Position (km)')
ax5[1].set_ylabel('Relative Velocity (km/s)')
ax5[2].set_ylabel('Deceleration (km/s$^2$)')
plt.show()

The above plot shows the deceleration values are inversely proportional to altitude.

Negative deceleration may be due to the orbit and resulting relative velocities w.r.t earth rotation.

In [159]:
%matplotlib notebook
#fig = plt.figure(dpi=100) #
fig = plt.figure(figsize=(6,6),dpi=135)
#fig = plt.figure(dpi=100)
#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
ax.set_xlim([-8000,8000])
ax.set_ylim([-8000,8000])
ax.set_zlim([-8000,8000])

###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)

#Show  Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation: Earth')
ax.scatter(0,0,0,color='k',alpha=0.8)

###########################################
#Plot First Full Orbit
clip1 = 0
clip2 = 7500
ax.plot(r1[clip1:clip2,0], r1[clip1:clip2,1], r1[clip1:clip2,2],'k',lw=1,zorder = 2.5,label='Initial Orbit')

###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'

i=0
j1 = r1_apo[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_ra,s=30) #show location of apoapsis
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_ra,label='$r_{Apoapsis}$') # show vector to apoapsis

#Show velocity vector at apoapsis
scale1 = 1000
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v)

#Show decel vector at apoapsis
scale2 = 1e25
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale2,a1[j1,1]*scale2,a1[j1,2]*scale2, color=c_a)

###########################################
#perigee
i=0
j1 = r1_peri[0][i]
ax.scatter(r1[j1,0],r1[j1,1],r1[j1,2], c=c_rp,s=30)#show location of perigee
ax.quiver(0,0,0,r1[j1,0],r1[j1,1],r1[j1,2], color=c_rp,label='$r_{Periapsis}$') # show vector to perigee

#Show velocity vector at perigee
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],v1rel[j1,0]*scale1,v1rel[j1,1]*scale1,v1rel[j1,2]*scale1, color=c_v,label='$\\vec v$')

#Show deceleration vector at perigee (Scales of deceleration are not equal)
scale3 = 1e20
ax.quiver(r1[j1,0],r1[j1,1],r1[j1,2],a1[j1,0]*scale3,a1[j1,1]*scale3,a1[j1,2]*scale3, color=c_a,label='$ \\vec a$ (not to scale)')

###########################################
#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')

#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])


plt.legend()
ax.set_axis_off()
  
plt.show()    

animate_rotate = True
if animate_rotate==True:
    for angle in range(0, 360):
        ax.view_init(0, angle)
        fig.canvas.draw()
        plt.pause(.001)

ax.view_init(-22,70)

The above figure is another representation of the inverse proportionality between deceleration and altitude.

DO AGAIN BUT ANIMATE THE TRAJECTORY

In [181]:
%matplotlib notebook
#fig = plt.figure(dpi=100) #
#fig = plt.figure(dpi=100)


#Very good
fig = plt.figure(figsize=(6,6),dpi=135)



#ax = fig.add_subplot(111, projection = '3d')
ax = Axes3D(fig)
axlim = 7000
ax.set_xlim([-axlim,axlim])
ax.set_ylim([-axlim,axlim])
ax.set_zlim([-axlim,axlim])
ax.view_init(-22,70)

###########################################
#Add Earth (Bottom Half)
ue1, ve1 = np.mgrid[0:2*np.pi:60j, 0:np.pi:60j]
xe=re*np.cos(ue1)*np.sin(ve1)
ye=re*np.sin(ue1)*np.sin(ve1)
ze=re*np.cos(ve1)
clip0 = int(xe.size / 10)
ax.plot_surface(xe, ye, ze, lw=0.1,alpha=0.2,color='deepskyblue',zorder = -10)
ax.plot_wireframe(xe, ye, ze, lw=0.1,alpha=0.99,color='deepskyblue',zorder = -10)

#Show  Earth Axis of rotation.
ax.plot([0,0],[0,0],[-re*1.2,re*1.2],color='k',lw=4,ls='-',alpha=0.2,label='Axis of Rotation')
ax.scatter(0,0,0,color='k',alpha=0.8)

###########################################
#Plot First Full Orbit
#clip1 = 0
#clip2 = 7500
#

###########################################
#apoapsis
c_v = 'mediumspringgreen'
c_a = 'hotpink'


#Label
ax.set_xlabel('x (km)')
ax.set_ylabel('y (km)')
ax.set_zlabel('z (km)')

#Set ticks
ax.set_xticks([-5000,0,5000])
ax.set_yticks([-5000,0,5000])
ax.set_zticks([-5000,0,5000])


#ax.set_axis_off()
def clear_lines():
    quiv_accel.remove()
    quiv_vel.remove()
    quiv_radius.remove()
    

plt.show()

LW_vec = 3


#For parameter output in plot
line = -20000
ypos_0 = 0.7 #-0.1
xpos_0 =0.02
    
    
animation_range = 1000
for j in range(animation_range):    

    
    i = j * 30   # * 80 for display
    
    #ORBIT
    ax.plot((r1[i,0],r1[i+1,0]), (r1[i,1],r1[i+1,1]), (r1[i,2],r1[i+1,2]),'k',lw=1,color='k',zorder = 2.5,label='Orbit')
    
    #RADIUS VECTOR
    quiv_radius = ax.quiver(0,0,0,r1[i,0],r1[i,1],r1[i,2], color=c_rp,linewidths=LW_vec, label='Vector to sat') 
    
    #VELOCITY VECTOR
    s1 = 1000
    quiv_vel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],v1rel[i,0]*s1,v1rel[i,1]*s1,v1rel[i,2]*s1,color=c_v,linewidths=LW_vec,label='Velocity')
    
    #DECEL VECTOR
    s2 = 1e20
    quiv_accel = ax.quiver(r1[i,0],r1[i,1],r1[i,2],a1[i,0]*s2,a1[i,1]*s2,a1[i,2]*s2, color=c_a, linewidths=LW_vec,label = 'Drag')
    
    
    
    text_below = '$altitude =$'+'{0:.0f} km \ndrag = {1:.1e} $km/s^2$'.format(r1_mag[i] - re,a1_mag[i]  )      
    text_parameter_1 = ax.text(xpos_0, ypos_0 + 0*line ,0,text_below,
        horizontalalignment='left',
        verticalalignment='bottom',
        size='large',
        bbox=dict(facecolor='gray', alpha=0.3),
        transform = ax.transAxes)
    
    

    
    #SHOW
    fig.canvas.draw()
    plt.pause(0.0001)

    if j == 0:
        plt.legend(loc='lower right',  borderaxespad=0.) #bbox_to_anchor=(1.05, 1),
    
    #Save frames for video
    frame_name = 'animation/frame_{}'.format(j)
    plt.savefig(frame_name)
    
    #REMOVE PREVIOUS
    if j < animation_range-1:
        quiv_accel.remove()
        quiv_vel.remove()
        quiv_radius.remove()
        text_parameter_1.remove()

        del quiv_accel
        del quiv_vel
        del quiv_radius
        del text_parameter_1

    #quiv_vel.set_alpha(0)
    #quiv_radius.set_alpha(0)
    #quiv_accel.set_alpha(0)
    #

###########################################
  
plt.show()